L’objectif de ce TP est de traiter le cas d’un mélange gaussien (à deux composantes pour le moment).

Simulation d’un mélange

On peut commencer par simuler un mélange pour travailler avec des données “faciles” :

rmix = function(N,p,moy1,moy2,t1,t2){
  b = rbinom(N,1,p)
  x = (1-b)*rnorm(N,moy1,1/sqrt(t1))+(b)*rnorm(N,moy2,1/sqrt(t2))
  return(x)
  }
hist(rmix(1000,0.6,0,5,1,1),100,main = "Simulated Mixture")

Modèle à variables cachés

On travaille maintenant avec le modèle suivant :

Xi|Zi,m1,m2,t1,t2i.i.dN(mZi,1tZi)Zi|pi.i.d1+B(p)pBeta(a0,b0)miN(μ0,σ02),i=1,2τiGamma(k0,λ0),i=1,2

On peut initialiser les données :

N = 1000
true = list(p = 0.6,m1 = 0,m2 = 4,tau1 =1,tau2 = 5)##pour la simulation
x = rmix(1000,true$p,true$m1,true$m2,true$tau1,true$tau2)
h = list(a =1,b = 1, mu = 0,sig = 1000,k = 1,lam = 0.001)
l = list(p  = 1/2,m = c(0,0),tau = c(1,1),z = 1+rbinom(length(x),1,1/2))
init = function(x){
  h = list(a =1,b = 1, mu = 0,sig = 1000,k = 1,lam = 0.001)
  l = list(p  = 1/2,m = c(0,0),tau = c(1,1),z = 1+rbinom(length(x),1,1/2))
  return(list(h = h,l = l))
}
logvraiss= function(l,x){
  return( sum( log( dnorm(x,l$m[l$z],1/sqrt(l$tau[l$z])    )   )       )         ) 
}
logvraiss(l,x)
[1] -5851.625
logprior = function(l,h){
  val = sum(log(dbinom(l$z-1,1,l$p )))+
    log(dbeta(l$p,h$a,h$b))+
    sum(log(dnorm(l$m,h$mu,h$sig)  ))+
    sum(log(dgamma(l$tau,h$k,h$lam)))
  return(val)
}
logprior(l,h)
[1] -722.6181
logpost = function(l,x,h){
    return(logvraiss(l,x)+logprior(l,h))
}
logpost(l,x,h)
[1] -6574.243
step_MH = function(l,x,h,sigprop = c(0.01,0.01,0.01),nrunz = 50){
  ##variance prop codées en dur pour le moment
  sig_prop_m = sigprop[1]
  inv_prop_tau = 1/sigprop[2]
  inv_prop_p = 1/sigprop[3]
  lnew = l
  N = length(x)
  # print(logpost(lnew,x,h))
  ##ici, c'est bien couteux de recalculer tout à chaque fois -> on devrait calculer que pour x[i] la vraisemblance !!!!
  ##proposition m
  i = sample(2,1)
  lnew$m[i] = lnew$m[i]+rnorm(1,0,sig_prop_m)
  if (log(runif(1))>(logpost(lnew,x,h)-logpost(l,x,h)) ){
    ##on rejette avec proba 1-alpha,  ici la proposition est symétrique
    lnew$m[i] = l$m[i]
  }
  # print(logpost(lnew,x,h))
  if (lnew$m[2]<lnew$m[1]){
    ##si les m se sont croisés, on échange tout : m, tau, p et z
    lnew$m = rev(lnew$m)
    lnew$tau = rev(lnew$tau)
    lnew$p = 1-lnew$p
    lnew$z = 3-lnew$z
  }
  # print(logpost(lnew,x,h))
  ##proposition tau
  i = sample(2,1)
  lnew$tau[i] = rgamma(1,inv_prop_tau*l$tau[i]**2,inv_prop_tau *l$tau[i])
  log_ratio_prop = log(dgamma(l$tau[i]     ,inv_prop_tau*lnew$tau[i]**2,inv_prop_tau *lnew$tau[i]))-log(dgamma(lnew$tau[i]     ,inv_prop_tau*l$tau[i]**2,inv_prop_tau *l$tau[i]))
  if (log(runif(1))>(log_ratio_prop+logpost(lnew,x,h)-logpost(l,x,h)) ){
    ##on rejette avec proba 1-alpha. ici on prend en compte le ratio
    lnew$tau[i] = l$tau[i]
  }
  ##proposition z
  # for (i in 1:N){
  # nrunz = 50
  samp = sample(N,nrunz)
  for (i in samp){
    # if (runif(1)<1/2){
      # i = sample(N,1)
      lnew$z[i] = 3-l$z[i]
      if (log(runif(1))>(logpost(lnew,x,h)-logpost(l,x,h)) ){
        ##on rejette avec proba 1-alpha. ici la proposition est symétrique
        lnew$z[i] = l$z[i]
      } else{
        l$z[i] = lnew$z[i]
      }
    # }
  }
  # print(logpost(lnew,x,h))
  ##proposition p
  lnew$p = rbeta(1,inv_prop_p*l$p,inv_prop_p*(1-l$p))
  log_ratio_prop = log(dbeta(l$p,inv_prop_p*lnew$p,inv_prop_p*(1-lnew$p)))-log(dbeta(lnew$p,inv_prop_p*l$p,inv_prop_p*(1-l$p)))
    if (log(runif(1))>(log_ratio_prop+logpost(lnew,x,h)-logpost(l,x,h)) ){
    ##on rejette avec proba 1-alpha. ici on prend en compte le ratio
    lnew$p = l$p
    }
  return(lnew)
}
lnew = step_MH(l,x,h)
run_MH = function(nMH,x,...){
  linit = init(x)
  l = linit$l
  h = linit$h
  traj = matrix(0,nMH,5)
  colnames(traj) = c("p","m0","m1","tau0","tau1")
  nb = rep(0,length(x))
  for (i in 1:nMH){
    nb = nb+(l$z-1)
    l = step_MH(l,x,h,...)
    traj[i,] = c(l$p,l$m,l$tau)
  }
  class = nb/nMH
  return(list(traj = traj ,class =class,l = l))
}
nMH = 100
ltraj  = run_MH(nMH,x)
plot(1:nMH,ltraj$traj[,2],"l")
lines(1:nMH,ltraj$traj[,3],col = "red")

Pour un plot plus joli, on reprend la fonction de la dernière fois

plot_run = function(x,mem,l,i = NULL,nbpoints = 10000){
    if (is.null(i)){i = dim(mem)[1]}
    par(mfrow = c(1,2))
    hist(x,50,freq= F,main = "Current estimated density")
    curve((1-l$p)*dnorm(x,l$m1,1/sqrt(l$tau1))+(l$p)*dnorm(x,l$m2,1/sqrt(l$tau2)),add = T,col = "red",lty = 2,lwd=  2)
    curve((1-l$p)*dnorm(x,l$m1,1/sqrt(l$tau1)),add = T,col = "blue",lty = 3,lwd = 2)
    curve((l$p)*dnorm(x,l$m2,1/sqrt(l$tau2)),add = T,col = "green",lty = 3,lwd = 2)
    z1 = sample(i,nbpoints,replace = T)
    z2 = sample(i,nbpoints,replace = T)
    plot(z1,mem[z1,2]+2*1/sqrt(mem[z1,4])*runif(nbpoints,-1,1),ylim = range(x),col = "cyan",cex = 0.1,main = "MH trajectory with 2-sigma bands",ylab = "Trajectories for m1,m2",xlab= "Time")
    points(z2,mem[z2,3]+2*1/sqrt(mem[z2,5])*runif(nbpoints,-1,1),ylim = range(x),col = "yellow",cex = 0.1)
    lines(1:i,mem[1:i,2],ylim = range(x),col = "blue")
    lines(1:i,mem[1:i,3],col = "green")
}
init = function(x){
  h = list(a =1,b = 1, mu = 0,sig = 1000,k = 1,lam = 0.001)
  l = list(p  = 1/2,m = c(0,0),tau = c(1,1),z = 1+rbinom(length(x),1,1/2))
  # l = list(p  = 1/2,m = c(-1,1),tau = c(1,1),z = 1+(x>mean(x)))
  return(list(h = h,l = l))
}
N = 1000
true = list(p = 0.6,m1 = 0,m2 = 2,tau1 =1,tau2 = 5)##pour la simulation
x = rmix(1000,true$p,true$m1,true$m2,true$tau1,true$tau2)
nMH = 2000
sigprop =  c(0.5,0.1,0.01)
nrunz = 50
ltraj  = run_MH(nMH,x,sigprop,nrunz)
l = ltraj$l
lend = list(p = l$p,m1  = l$m[1],m2 = l$m[2],tau1 = l$tau[1],tau2 = l$tau[2])
plot_run(x,ltraj$traj,lend)

On peut alors voir les paramètres finaux :

cat("estimés : ",unlist(lend),"\n")
estimés :  0.5321751 0.4869726 1.986869 0.7883785 5.634391 
cat("théorique :",unlist(true),"\n")
théorique : 0.6 0 2 1 5 

On peut aussi observer la propotion du temps que chaque observation a passé avec Zi=1 :

s = sort(x,index.return = T)
plot(s$x,ltraj$class[s$ix],"l",main  = "Classification probability",xlab = "x",ylab = "Estimated P(Z(x)= 1)")

Et si on veut faire un petit film, il n’y a qu’à modifier les fonctions

library(animation)
plot_run = function(x,mem,class,l,i = NULL,nbpoints = 10000,MHtot = NULL){
    if (is.null(i)){i = dim(mem)[1]}
    if (is.null(MHtot)){MHtot = dim(mem)[1]}
    # par(mfrow = c(1,3))
    hist(x,50,freq= F,main = "Current estimated density")
    curve((1-l$p)*dnorm(x,l$m1,1/sqrt(l$tau1))+(l$p)*dnorm(x,l$m2,1/sqrt(l$tau2)),add = T,col = "red",lty = 2,lwd=  3)
    curve((1-l$p)*dnorm(x,l$m1,1/sqrt(l$tau1)),add = T,col = "blue",lty = 3,lwd = 3)
    curve((l$p)*dnorm(x,l$m2,1/sqrt(l$tau2)),add = T,col = "green",lty = 3,lwd = 3)
    z1 = sample(i,nbpoints,replace = T)
    z2 = sample(i,nbpoints,replace = T)
    plot(z1,mem[z1,2]+2*1/sqrt(mem[z1,4])*runif(nbpoints,-1,1),ylim = range(x),xlim = c(1,MHtot),col = "cyan",cex = 0.1,main = "MH trajectory with 2-sigma bands")
    points(z2,mem[z2,3]+2*1/sqrt(mem[z2,5])*runif(nbpoints,-1,1),ylim = range(x),col = "yellow",cex = 0.1)
    lines(1:i,mem[1:i,2],ylim = range(x),col = "blue")
    lines(1:i,mem[1:i,3],col = "green")
    s = sort(x,index.return = T)
    plot(s$x,class[s$ix]/i,"l")
}
run_MH = function(nMH,x,pas_anim,...){
  linit = init(x)
  l = linit$l
  h = linit$h
  traj = matrix(0,nMH,5)
  colnames(traj) = c("p","m0","m1","tau0","tau1")
  nb = rep(0,length(x))
  traj[1,] = c(l$p,l$m,l$tau)
  for (i in 2:nMH){
    nb = nb+(l$z-1)
    l = step_MH(l,x,h,...)
    lnow = list(p = l$p,m1  = l$m[1],m2 = l$m[2],tau1 = l$tau[1],tau2 = l$tau[2])
    if ((i%%pas_anim)==2){
      plot_run(x,traj[1:i,],class = nb,l = lnow,i=i,MHtot =nMH )
      ani.record()
    }
    traj[i,] = c(l$p,l$m,l$tau)
  }
  class = nb/nMH
  return(list(traj = traj ,class =class,l = l))
}
# 
# N = 100
# true = list(p = 0.6,m1 = 0,m2 = 2,tau1 =1,tau2 = 5)##pour la simulation
# x = rmix(1000,true$p,true$m1,true$m2,true$tau1,true$tau2)
# sigprop =  0.1*c(0.5,0.1,0.01)
# nrunz = 20
# nMH = 500
# pas_anim = 10
# par(mfrow = c(1,3))
# ani.options(interval = 0.05, nmax = floor(nMH/pas_anim))
# ani.record(reset = TRUE)
# ltraj  = run_MH(nMH,x,pas_anim= pas_anim,sigprop,nrunz)
# saveGIF(ani.replay(),movie.name = "animation_MH.gif")

Commentaires et extensions

On peut aussi observer la chose suivante :

π(m1|X,Z,τ,p)π(m1)i=1N(exp(τ1(xim1)22))11zi=1exp((m1μ0)22σ02i=1N11zi=1τ1(xim1)22)

Par propriété de conjugaison de loi, on voit alors que

m1|X,Z,τ,pN(μ0σ02+τ1i=1N11zi=1xi1σ02+τ1i=1N11zi=1,11σ02+τ1i=1N11zi=1)

Et de même pour tous les autres paramètres. Cela incite à reprendre tout notre code pour calculer seulement cela à chaque étape, sans devoir recalculer la vraisemblance.

Voilà un bon exercice pour votre projet !

Et avec les packages

setwd("/home/espinasse/Bureau/TAF/Seafile/Cours/stat bayesiennes/tp R 2020/old_site")
library(rjags)
Le chargement a nécessité le package : coda
Linked to JAGS 4.3.0
Loaded modules: basemod,bugs
N = 1000
true = list(p = 0.6,m1 = 0,m2 = 2,tau1 =1,tau2 = 5)##pour la simulation
u = rmix(10000,true$p,true$m1,true$m2,true$tau1,true$tau2)
m <- jags.model("model_mixture.bug", data=list(x = u, N=N))
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 1000
   Unobserved stochastic nodes: 1005
   Total graph size: 14010

Initializing model
f = coda.samples(m,c("m","tau","p"),n.iter = 1000)

  |                                                        
  |                                                  |   0%
  |                                                        
  |*                                                 |   2%
  |                                                        
  |**                                                |   4%
  |                                                        
  |***                                               |   6%
  |                                                        
  |****                                              |   8%
  |                                                        
  |*****                                             |  10%
  |                                                        
  |******                                            |  12%
  |                                                        
  |*******                                           |  14%
  |                                                        
  |********                                          |  16%
  |                                                        
  |*********                                         |  18%
  |                                                        
  |**********                                        |  20%
  |                                                        
  |***********                                       |  22%
  |                                                        
  |************                                      |  24%
  |                                                        
  |*************                                     |  26%
  |                                                        
  |**************                                    |  28%
  |                                                        
  |***************                                   |  30%
  |                                                        
  |****************                                  |  32%
  |                                                        
  |*****************                                 |  34%
  |                                                        
  |******************                                |  36%
  |                                                        
  |*******************                               |  38%
  |                                                        
  |********************                              |  40%
  |                                                        
  |*********************                             |  42%
  |                                                        
  |**********************                            |  44%
  |                                                        
  |***********************                           |  46%
  |                                                        
  |************************                          |  48%
  |                                                        
  |*************************                         |  50%
  |                                                        
  |**************************                        |  52%
  |                                                        
  |***************************                       |  54%
  |                                                        
  |****************************                      |  56%
  |                                                        
  |*****************************                     |  58%
  |                                                        
  |******************************                    |  60%
  |                                                        
  |*******************************                   |  62%
  |                                                        
  |********************************                  |  64%
  |                                                        
  |*********************************                 |  66%
  |                                                        
  |**********************************                |  68%
  |                                                        
  |***********************************               |  70%
  |                                                        
  |************************************              |  72%
  |                                                        
  |*************************************             |  74%
  |                                                        
  |**************************************            |  76%
  |                                                        
  |***************************************           |  78%
  |                                                        
  |****************************************          |  80%
  |                                                        
  |*****************************************         |  82%
  |                                                        
  |******************************************        |  84%
  |                                                        
  |*******************************************       |  86%
  |                                                        
  |********************************************      |  88%
  |                                                        
  |*********************************************     |  90%
  |                                                        
  |**********************************************    |  92%
  |                                                        
  |***********************************************   |  94%
  |                                                        
  |************************************************  |  96%
  |                                                        
  |************************************************* |  98%
  |                                                        
  |**************************************************| 100%
traj = f[[1]]
summary(f)

Iterations = 1:1000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 1000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

           Mean      SD  Naive SE Time-series SE
m[1]   -0.02481 0.14448 0.0045689       0.018585
m[2]    1.95753 0.05871 0.0018567       0.005077
p       0.62579 0.02915 0.0009219       0.003282
tau[1]  0.92491 0.13191 0.0041714       0.015436
tau[2]  5.39071 0.60625 0.0191714       0.064464

2. Quantiles for each variable:

          2.5%     25%      50%     75%  97.5%
m[1]   -0.2420 -0.1075 -0.03618 0.04833 0.1887
m[2]    1.9100  1.9480  1.96294 1.97752 2.0037
p       0.5721  0.6085  0.62655 0.64461 0.6777
tau[1]  0.6953  0.8372  0.91221 1.01074 1.2072
tau[2]  4.5220  5.1226  5.40535 5.71578 6.3382
par(mfrow = c(5,2))
# for (i in 1:5){
#   plot(1:nrow(traj),traj[,i],"l",main = colnames(traj)[i])
# }
plot(f)

LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKTCdvYmplY3RpZiBkZSBjZSBUUCBlc3QgZGUgdHJhaXRlciBsZSBjYXMgZCd1biBtw6lsYW5nZSBnYXVzc2llbiAow6AgZGV1eCBjb21wb3NhbnRlcyBwb3VyIGxlIG1vbWVudCkuIAoKCiMgU2ltdWxhdGlvbiBkJ3VuIG3DqWxhbmdlCk9uIHBldXQgY29tbWVuY2VyIHBhciBzaW11bGVyIHVuIG3DqWxhbmdlIHBvdXIgdHJhdmFpbGxlciBhdmVjIGRlcyBkb25uw6llcyAiZmFjaWxlcyIgOgoKYGBge3J9CnJtaXggPSBmdW5jdGlvbihOLHAsbW95MSxtb3kyLHQxLHQyKXsKICBiID0gcmJpbm9tKE4sMSxwKQogIHggPSAoMS1iKSpybm9ybShOLG1veTEsMS9zcXJ0KHQxKSkrKGIpKnJub3JtKE4sbW95MiwxL3NxcnQodDIpKQogIHJldHVybih4KQogIH0KaGlzdChybWl4KDEwMDAsMC42LDAsNSwxLDEpLDEwMCxtYWluID0gIlNpbXVsYXRlZCBNaXh0dXJlIikKYGBgCgoKIyBNb2TDqGxlIMOgIHZhcmlhYmxlcyBjYWNow6lzIAoKCk9uIHRyYXZhaWxsZSBtYWludGVuYW50IGF2ZWMgbGUgbW9kw6hsZSBzdWl2YW50IDoKCiQkXGJlZ2lue2FsaWduKn0KWF9pfFpfaSxtXzEsbV8yLHRfMSx0XzImIFxzaW1fe2kuaS5kfSBcbWF0aGNhbHtOfVxCaWcobV97Wl9pfSxcZnJhY3sxfXt0X3taX2l9fVxCaWcpClxcWl9pfHAgJlxzaW1fe2kuaS5kfSAxK1xtYXRoY2Fse0J9KHApClxcIHAgJlxzaW0gQmV0YShhXzAsYl8wKQpcXCBtX2kgJlxzaW0gXG1hdGhjYWx7Tn1cQmlnKFxtdV8wLFxzaWdtYV8wXjJcQmlnKSwgaSA9IDEsIDIKXFwgXHRhdV9pICZcc2ltIEdhbW1hKGtfMCxcbGFtYmRhXzApLCBpID0gMSwyClxlbmR7YWxpZ24qfSQkCgoKT24gcGV1dCBpbml0aWFsaXNlciBsZXMgZG9ubsOpZXMgOgpgYGB7cn0KTiA9IDEwMDAKdHJ1ZSA9IGxpc3QocCA9IDAuNixtMSA9IDAsbTIgPSA0LHRhdTEgPTEsdGF1MiA9IDUpIyNwb3VyIGxhIHNpbXVsYXRpb24KeCA9IHJtaXgoMTAwMCx0cnVlJHAsdHJ1ZSRtMSx0cnVlJG0yLHRydWUkdGF1MSx0cnVlJHRhdTIpCmggPSBsaXN0KGEgPTEsYiA9IDEsIG11ID0gMCxzaWcgPSAxMDAwLGsgPSAxLGxhbSA9IDAuMDAxKQpsID0gbGlzdChwICA9IDEvMixtID0gYygwLDApLHRhdSA9IGMoMSwxKSx6ID0gMStyYmlub20obGVuZ3RoKHgpLDEsMS8yKSkKYGBgCgoKYGBge3J9Cgppbml0ID0gZnVuY3Rpb24oeCl7CiAgaCA9IGxpc3QoYSA9MSxiID0gMSwgbXUgPSAwLHNpZyA9IDEwMDAsayA9IDEsbGFtID0gMC4wMDEpCiAgbCA9IGxpc3QocCAgPSAxLzIsbSA9IGMoMCwwKSx0YXUgPSBjKDEsMSkseiA9IDErcmJpbm9tKGxlbmd0aCh4KSwxLDEvMikpCiAgcmV0dXJuKGxpc3QoaCA9IGgsbCA9IGwpKQp9Cgpsb2d2cmFpc3M9IGZ1bmN0aW9uKGwseCl7CiAgcmV0dXJuKCBzdW0oIGxvZyggZG5vcm0oeCxsJG1bbCR6XSwxL3NxcnQobCR0YXVbbCR6XSkgICAgKSAgICkgICAgICAgKSAgICAgICAgICkgCn0KCmxvZ3ZyYWlzcyhsLHgpCgpsb2dwcmlvciA9IGZ1bmN0aW9uKGwsaCl7CiAgdmFsID0gc3VtKGxvZyhkYmlub20obCR6LTEsMSxsJHAgKSkpKwogICAgbG9nKGRiZXRhKGwkcCxoJGEsaCRiKSkrCiAgICBzdW0obG9nKGRub3JtKGwkbSxoJG11LGgkc2lnKSAgKSkrCiAgICBzdW0obG9nKGRnYW1tYShsJHRhdSxoJGssaCRsYW0pKSkKICByZXR1cm4odmFsKQp9Cgpsb2dwcmlvcihsLGgpCgpsb2dwb3N0ID0gZnVuY3Rpb24obCx4LGgpewogICAgcmV0dXJuKGxvZ3ZyYWlzcyhsLHgpK2xvZ3ByaW9yKGwsaCkpCn0KCmxvZ3Bvc3QobCx4LGgpCmBgYAoKCgoKYGBge3J9CnN0ZXBfTUggPSBmdW5jdGlvbihsLHgsaCxzaWdwcm9wID0gYygwLjAxLDAuMDEsMC4wMSksbnJ1bnogPSA1MCl7CiAgIyN2YXJpYW5jZSBwcm9wIGNvZMOpZXMgZW4gZHVyIHBvdXIgbGUgbW9tZW50CiAgc2lnX3Byb3BfbSA9IHNpZ3Byb3BbMV0KICBpbnZfcHJvcF90YXUgPSAxL3NpZ3Byb3BbMl0KICBpbnZfcHJvcF9wID0gMS9zaWdwcm9wWzNdCiAgbG5ldyA9IGwKICBOID0gbGVuZ3RoKHgpCiAgIyBwcmludChsb2dwb3N0KGxuZXcseCxoKSkKICAjI2ljaSwgYydlc3QgYmllbiBjb3V0ZXV4IGRlIHJlY2FsY3VsZXIgdG91dCDDoCBjaGFxdWUgZm9pcyAtPiBvbiBkZXZyYWl0IGNhbGN1bGVyIHF1ZSBwb3VyIHhbaV0gbGEgdnJhaXNlbWJsYW5jZSAhISEhCiAgIyNwcm9wb3NpdGlvbiBtCiAgaSA9IHNhbXBsZSgyLDEpCiAgbG5ldyRtW2ldID0gbG5ldyRtW2ldK3Jub3JtKDEsMCxzaWdfcHJvcF9tKQogIGlmIChsb2cocnVuaWYoMSkpPihsb2dwb3N0KGxuZXcseCxoKS1sb2dwb3N0KGwseCxoKSkgKXsKICAgICMjb24gcmVqZXR0ZSBhdmVjIHByb2JhIDEtYWxwaGEsICBpY2kgbGEgcHJvcG9zaXRpb24gZXN0IHN5bcOpdHJpcXVlCiAgICBsbmV3JG1baV0gPSBsJG1baV0KICB9CiAgIyBwcmludChsb2dwb3N0KGxuZXcseCxoKSkKICBpZiAobG5ldyRtWzJdPGxuZXckbVsxXSl7CiAgICAjI3NpIGxlcyBtIHNlIHNvbnQgY3JvaXPDqXMsIG9uIMOpY2hhbmdlIHRvdXQgOiBtLCB0YXUsIHAgZXQgegogICAgbG5ldyRtID0gcmV2KGxuZXckbSkKICAgIGxuZXckdGF1ID0gcmV2KGxuZXckdGF1KQogICAgbG5ldyRwID0gMS1sbmV3JHAKICAgIGxuZXckeiA9IDMtbG5ldyR6CiAgfQogICMgcHJpbnQobG9ncG9zdChsbmV3LHgsaCkpCiAgIyNwcm9wb3NpdGlvbiB0YXUKICBpID0gc2FtcGxlKDIsMSkKICBsbmV3JHRhdVtpXSA9IHJnYW1tYSgxLGludl9wcm9wX3RhdSpsJHRhdVtpXSoqMixpbnZfcHJvcF90YXUgKmwkdGF1W2ldKQogIGxvZ19yYXRpb19wcm9wID0gbG9nKGRnYW1tYShsJHRhdVtpXSAgICAgLGludl9wcm9wX3RhdSpsbmV3JHRhdVtpXSoqMixpbnZfcHJvcF90YXUgKmxuZXckdGF1W2ldKSktbG9nKGRnYW1tYShsbmV3JHRhdVtpXSAgICAgLGludl9wcm9wX3RhdSpsJHRhdVtpXSoqMixpbnZfcHJvcF90YXUgKmwkdGF1W2ldKSkKICBpZiAobG9nKHJ1bmlmKDEpKT4obG9nX3JhdGlvX3Byb3ArbG9ncG9zdChsbmV3LHgsaCktbG9ncG9zdChsLHgsaCkpICl7CiAgICAjI29uIHJlamV0dGUgYXZlYyBwcm9iYSAxLWFscGhhLiBpY2kgb24gcHJlbmQgZW4gY29tcHRlIGxlIHJhdGlvCiAgICBsbmV3JHRhdVtpXSA9IGwkdGF1W2ldCiAgfQogICMjcHJvcG9zaXRpb24gegogICMgZm9yIChpIGluIDE6Til7CiAgIyBucnVueiA9IDUwCiAgc2FtcCA9IHNhbXBsZShOLG5ydW56KQogIGZvciAoaSBpbiBzYW1wKXsKICAgICMgaWYgKHJ1bmlmKDEpPDEvMil7CiAgICAgICMgaSA9IHNhbXBsZShOLDEpCiAgICAgIGxuZXckeltpXSA9IDMtbCR6W2ldCiAgICAgIGlmIChsb2cocnVuaWYoMSkpPihsb2dwb3N0KGxuZXcseCxoKS1sb2dwb3N0KGwseCxoKSkgKXsKICAgICAgICAjI29uIHJlamV0dGUgYXZlYyBwcm9iYSAxLWFscGhhLiBpY2kgbGEgcHJvcG9zaXRpb24gZXN0IHN5bcOpdHJpcXVlCiAgICAgICAgbG5ldyR6W2ldID0gbCR6W2ldCiAgICAgIH0gZWxzZXsKICAgICAgICBsJHpbaV0gPSBsbmV3JHpbaV0KICAgICAgfQogICAgIyB9CiAgfQogICMgcHJpbnQobG9ncG9zdChsbmV3LHgsaCkpCiAgIyNwcm9wb3NpdGlvbiBwCiAgbG5ldyRwID0gcmJldGEoMSxpbnZfcHJvcF9wKmwkcCxpbnZfcHJvcF9wKigxLWwkcCkpCiAgbG9nX3JhdGlvX3Byb3AgPSBsb2coZGJldGEobCRwLGludl9wcm9wX3AqbG5ldyRwLGludl9wcm9wX3AqKDEtbG5ldyRwKSkpLWxvZyhkYmV0YShsbmV3JHAsaW52X3Byb3BfcCpsJHAsaW52X3Byb3BfcCooMS1sJHApKSkKICAgIGlmIChsb2cocnVuaWYoMSkpPihsb2dfcmF0aW9fcHJvcCtsb2dwb3N0KGxuZXcseCxoKS1sb2dwb3N0KGwseCxoKSkgKXsKICAgICMjb24gcmVqZXR0ZSBhdmVjIHByb2JhIDEtYWxwaGEuIGljaSBvbiBwcmVuZCBlbiBjb21wdGUgbGUgcmF0aW8KICAgIGxuZXckcCA9IGwkcAogICAgfQogIHJldHVybihsbmV3KQp9CmxuZXcgPSBzdGVwX01IKGwseCxoKQoKYGBgCgoKCgoKYGBge3J9CnJ1bl9NSCA9IGZ1bmN0aW9uKG5NSCx4LC4uLil7CiAgbGluaXQgPSBpbml0KHgpCiAgbCA9IGxpbml0JGwKICBoID0gbGluaXQkaAogIHRyYWogPSBtYXRyaXgoMCxuTUgsNSkKICBjb2xuYW1lcyh0cmFqKSA9IGMoInAiLCJtMCIsIm0xIiwidGF1MCIsInRhdTEiKQogIG5iID0gcmVwKDAsbGVuZ3RoKHgpKQogIGZvciAoaSBpbiAxOm5NSCl7CiAgICBuYiA9IG5iKyhsJHotMSkKICAgIGwgPSBzdGVwX01IKGwseCxoLC4uLikKICAgIHRyYWpbaSxdID0gYyhsJHAsbCRtLGwkdGF1KQogIH0KICBjbGFzcyA9IG5iL25NSAogIHJldHVybihsaXN0KHRyYWogPSB0cmFqICxjbGFzcyA9Y2xhc3MsbCA9IGwpKQp9Cm5NSCA9IDEwMApsdHJhaiAgPSBydW5fTUgobk1ILHgpCnBsb3QoMTpuTUgsbHRyYWokdHJhalssMl0sImwiKQpsaW5lcygxOm5NSCxsdHJhaiR0cmFqWywzXSxjb2wgPSAicmVkIikKYGBgCgpQb3VyIHVuIHBsb3QgcGx1cyBqb2xpLCBvbiByZXByZW5kIGxhIGZvbmN0aW9uIGRlIGxhIGRlcm5pw6hyZSBmb2lzCmBgYHtyfQpwbG90X3J1biA9IGZ1bmN0aW9uKHgsbWVtLGwsaSA9IE5VTEwsbmJwb2ludHMgPSAxMDAwMCl7CiAgICBpZiAoaXMubnVsbChpKSl7aSA9IGRpbShtZW0pWzFdfQogICAgcGFyKG1mcm93ID0gYygxLDIpKQogICAgaGlzdCh4LDUwLGZyZXE9IEYsbWFpbiA9ICJDdXJyZW50IGVzdGltYXRlZCBkZW5zaXR5IikKICAgIGN1cnZlKCgxLWwkcCkqZG5vcm0oeCxsJG0xLDEvc3FydChsJHRhdTEpKSsobCRwKSpkbm9ybSh4LGwkbTIsMS9zcXJ0KGwkdGF1MikpLGFkZCA9IFQsY29sID0gInJlZCIsbHR5ID0gMixsd2Q9ICAyKQogICAgY3VydmUoKDEtbCRwKSpkbm9ybSh4LGwkbTEsMS9zcXJ0KGwkdGF1MSkpLGFkZCA9IFQsY29sID0gImJsdWUiLGx0eSA9IDMsbHdkID0gMikKICAgIGN1cnZlKChsJHApKmRub3JtKHgsbCRtMiwxL3NxcnQobCR0YXUyKSksYWRkID0gVCxjb2wgPSAiZ3JlZW4iLGx0eSA9IDMsbHdkID0gMikKICAgIHoxID0gc2FtcGxlKGksbmJwb2ludHMscmVwbGFjZSA9IFQpCiAgICB6MiA9IHNhbXBsZShpLG5icG9pbnRzLHJlcGxhY2UgPSBUKQogICAgcGxvdCh6MSxtZW1bejEsMl0rMioxL3NxcnQobWVtW3oxLDRdKSpydW5pZihuYnBvaW50cywtMSwxKSx5bGltID0gcmFuZ2UoeCksY29sID0gImN5YW4iLGNleCA9IDAuMSxtYWluID0gIk1IIHRyYWplY3Rvcnkgd2l0aCAyLXNpZ21hIGJhbmRzIix5bGFiID0gIlRyYWplY3RvcmllcyBmb3IgbTEsbTIiLHhsYWI9ICJUaW1lIikKICAgIHBvaW50cyh6MixtZW1bejIsM10rMioxL3NxcnQobWVtW3oyLDVdKSpydW5pZihuYnBvaW50cywtMSwxKSx5bGltID0gcmFuZ2UoeCksY29sID0gInllbGxvdyIsY2V4ID0gMC4xKQogICAgbGluZXMoMTppLG1lbVsxOmksMl0seWxpbSA9IHJhbmdlKHgpLGNvbCA9ICJibHVlIikKICAgIGxpbmVzKDE6aSxtZW1bMTppLDNdLGNvbCA9ICJncmVlbiIpCn0KCmluaXQgPSBmdW5jdGlvbih4KXsKICBoID0gbGlzdChhID0xLGIgPSAxLCBtdSA9IDAsc2lnID0gMTAwMCxrID0gMSxsYW0gPSAwLjAwMSkKICBsID0gbGlzdChwICA9IDEvMixtID0gYygwLDApLHRhdSA9IGMoMSwxKSx6ID0gMStyYmlub20obGVuZ3RoKHgpLDEsMS8yKSkKICAjIGwgPSBsaXN0KHAgID0gMS8yLG0gPSBjKC0xLDEpLHRhdSA9IGMoMSwxKSx6ID0gMSsoeD5tZWFuKHgpKSkKICByZXR1cm4obGlzdChoID0gaCxsID0gbCkpCn0KCk4gPSAxMDAwCnRydWUgPSBsaXN0KHAgPSAwLjYsbTEgPSAwLG0yID0gMix0YXUxID0xLHRhdTIgPSA1KSMjcG91ciBsYSBzaW11bGF0aW9uCnggPSBybWl4KDEwMDAsdHJ1ZSRwLHRydWUkbTEsdHJ1ZSRtMix0cnVlJHRhdTEsdHJ1ZSR0YXUyKQpuTUggPSAyMDAwCnNpZ3Byb3AgPSAgYygwLjUsMC4xLDAuMDEpCm5ydW56ID0gNTAKbHRyYWogID0gcnVuX01IKG5NSCx4LHNpZ3Byb3AsbnJ1bnopCmwgPSBsdHJhaiRsCmxlbmQgPSBsaXN0KHAgPSBsJHAsbTEgID0gbCRtWzFdLG0yID0gbCRtWzJdLHRhdTEgPSBsJHRhdVsxXSx0YXUyID0gbCR0YXVbMl0pCnBsb3RfcnVuKHgsbHRyYWokdHJhaixsZW5kKQoKYGBgCgpPbiBwZXV0IGFsb3JzIHZvaXIgbGVzIHBhcmFtw6h0cmVzIGZpbmF1eCA6CgpgYGB7cn0KY2F0KCJlc3RpbcOpcyA6ICIsdW5saXN0KGxlbmQpLCJcbiIpCmNhdCgidGjDqW9yaXF1ZSA6Iix1bmxpc3QodHJ1ZSksIlxuIikKYGBgCgpPbiBwZXV0IGF1c3NpIG9ic2VydmVyIGxhIHByb3BvdGlvbiBkdSB0ZW1wcyBxdWUgY2hhcXVlIG9ic2VydmF0aW9uIGEgcGFzc8OpIGF2ZWMgJFpfaT0xJCA6CgpgYGB7cn0KcyA9IHNvcnQoeCxpbmRleC5yZXR1cm4gPSBUKQpwbG90KHMkeCxsdHJhaiRjbGFzc1tzJGl4XSwibCIsbWFpbiAgPSAiQ2xhc3NpZmljYXRpb24gcHJvYmFiaWxpdHkiLHhsYWIgPSAieCIseWxhYiA9ICJFc3RpbWF0ZWQgUChaKHgpPSAxKSIpCmBgYAoKRXQgc2kgb24gdmV1dCBmYWlyZSB1biBwZXRpdCBmaWxtLCBpbCBuJ3kgYSBxdSfDoCBtb2RpZmllciBsZXMgZm9uY3Rpb25zCmBgYHtyfQpsaWJyYXJ5KGFuaW1hdGlvbikKCnBsb3RfcnVuID0gZnVuY3Rpb24oeCxtZW0sY2xhc3MsbCxpID0gTlVMTCxuYnBvaW50cyA9IDEwMDAwLE1IdG90ID0gTlVMTCl7CiAgICBpZiAoaXMubnVsbChpKSl7aSA9IGRpbShtZW0pWzFdfQogICAgaWYgKGlzLm51bGwoTUh0b3QpKXtNSHRvdCA9IGRpbShtZW0pWzFdfQogICAgIyBwYXIobWZyb3cgPSBjKDEsMykpCiAgICBoaXN0KHgsNTAsZnJlcT0gRixtYWluID0gIkN1cnJlbnQgZXN0aW1hdGVkIGRlbnNpdHkiKQogICAgY3VydmUoKDEtbCRwKSpkbm9ybSh4LGwkbTEsMS9zcXJ0KGwkdGF1MSkpKyhsJHApKmRub3JtKHgsbCRtMiwxL3NxcnQobCR0YXUyKSksYWRkID0gVCxjb2wgPSAicmVkIixsdHkgPSAyLGx3ZD0gIDMpCiAgICBjdXJ2ZSgoMS1sJHApKmRub3JtKHgsbCRtMSwxL3NxcnQobCR0YXUxKSksYWRkID0gVCxjb2wgPSAiYmx1ZSIsbHR5ID0gMyxsd2QgPSAzKQogICAgY3VydmUoKGwkcCkqZG5vcm0oeCxsJG0yLDEvc3FydChsJHRhdTIpKSxhZGQgPSBULGNvbCA9ICJncmVlbiIsbHR5ID0gMyxsd2QgPSAzKQogICAgejEgPSBzYW1wbGUoaSxuYnBvaW50cyxyZXBsYWNlID0gVCkKICAgIHoyID0gc2FtcGxlKGksbmJwb2ludHMscmVwbGFjZSA9IFQpCiAgICBwbG90KHoxLG1lbVt6MSwyXSsyKjEvc3FydChtZW1bejEsNF0pKnJ1bmlmKG5icG9pbnRzLC0xLDEpLHlsaW0gPSByYW5nZSh4KSx4bGltID0gYygxLE1IdG90KSxjb2wgPSAiY3lhbiIsY2V4ID0gMC4xLG1haW4gPSAiTUggdHJhamVjdG9yeSB3aXRoIDItc2lnbWEgYmFuZHMiKQogICAgcG9pbnRzKHoyLG1lbVt6MiwzXSsyKjEvc3FydChtZW1bejIsNV0pKnJ1bmlmKG5icG9pbnRzLC0xLDEpLHlsaW0gPSByYW5nZSh4KSxjb2wgPSAieWVsbG93IixjZXggPSAwLjEpCiAgICBsaW5lcygxOmksbWVtWzE6aSwyXSx5bGltID0gcmFuZ2UoeCksY29sID0gImJsdWUiKQogICAgbGluZXMoMTppLG1lbVsxOmksM10sY29sID0gImdyZWVuIikKICAgIHMgPSBzb3J0KHgsaW5kZXgucmV0dXJuID0gVCkKICAgIHBsb3QocyR4LGNsYXNzW3MkaXhdL2ksImwiKQp9CgoKcnVuX01IID0gZnVuY3Rpb24obk1ILHgscGFzX2FuaW0sLi4uKXsKICBsaW5pdCA9IGluaXQoeCkKICBsID0gbGluaXQkbAogIGggPSBsaW5pdCRoCiAgdHJhaiA9IG1hdHJpeCgwLG5NSCw1KQogIGNvbG5hbWVzKHRyYWopID0gYygicCIsIm0wIiwibTEiLCJ0YXUwIiwidGF1MSIpCiAgbmIgPSByZXAoMCxsZW5ndGgoeCkpCiAgdHJhalsxLF0gPSBjKGwkcCxsJG0sbCR0YXUpCiAgZm9yIChpIGluIDI6bk1IKXsKICAgIG5iID0gbmIrKGwkei0xKQogICAgbCA9IHN0ZXBfTUgobCx4LGgsLi4uKQogICAgbG5vdyA9IGxpc3QocCA9IGwkcCxtMSAgPSBsJG1bMV0sbTIgPSBsJG1bMl0sdGF1MSA9IGwkdGF1WzFdLHRhdTIgPSBsJHRhdVsyXSkKICAgIGlmICgoaSUlcGFzX2FuaW0pPT0yKXsKICAgICAgcGxvdF9ydW4oeCx0cmFqWzE6aSxdLGNsYXNzID0gbmIsbCA9IGxub3csaT1pLE1IdG90ID1uTUggKQogICAgICBhbmkucmVjb3JkKCkKICAgIH0KICAgIHRyYWpbaSxdID0gYyhsJHAsbCRtLGwkdGF1KQogIH0KICBjbGFzcyA9IG5iL25NSAogIHJldHVybihsaXN0KHRyYWogPSB0cmFqICxjbGFzcyA9Y2xhc3MsbCA9IGwpKQp9CiMgCiMgTiA9IDEwMAojIHRydWUgPSBsaXN0KHAgPSAwLjYsbTEgPSAwLG0yID0gMix0YXUxID0xLHRhdTIgPSA1KSMjcG91ciBsYSBzaW11bGF0aW9uCiMgeCA9IHJtaXgoMTAwMCx0cnVlJHAsdHJ1ZSRtMSx0cnVlJG0yLHRydWUkdGF1MSx0cnVlJHRhdTIpCiMgc2lncHJvcCA9ICAwLjEqYygwLjUsMC4xLDAuMDEpCiMgbnJ1bnogPSAyMAojIG5NSCA9IDUwMAojIHBhc19hbmltID0gMTAKIyBwYXIobWZyb3cgPSBjKDEsMykpCiMgYW5pLm9wdGlvbnMoaW50ZXJ2YWwgPSAwLjA1LCBubWF4ID0gZmxvb3Iobk1IL3Bhc19hbmltKSkKIyBhbmkucmVjb3JkKHJlc2V0ID0gVFJVRSkKIyBsdHJhaiAgPSBydW5fTUgobk1ILHgscGFzX2FuaW09IHBhc19hbmltLHNpZ3Byb3AsbnJ1bnopCgoKIyBzYXZlR0lGKGFuaS5yZXBsYXkoKSxtb3ZpZS5uYW1lID0gImFuaW1hdGlvbl9NSC5naWYiKQoKYGBgCgo8aW1nIHNyYz0iYW5pbWF0aW9uX01ILmdpZiI+CgoKCiMjIENvbW1lbnRhaXJlcyBldCBleHRlbnNpb25zCgoKT24gcGV1dCBhdXNzaSBvYnNlcnZlciBsYSBjaG9zZSBzdWl2YW50ZSA6CgokJFxiZWdpbnthbGlnbip9ClxwaShtXzF8WCxaLFx0YXUscCkgJlxwcm9wdG8gXHBpKG1fMSlccHJvZF97aSA9IDF9Xk4gXEJpZ2coXGV4cCBcYmlnZygtIFxmcmFje1x0YXVfMSh4X2ktbV8xKV4yfXsyfVxiaWdnKSBcQmlnZyleezFcIVwhMV97el9pID0gMX19ClxcICYgXHByb3B0byBcZXhwXEJpZ2coLVxmcmFjeyhtXzEtXG11XzApXjJ9ezJcc2lnbWFfMF4yfS0gXHN1bV97aSA9IDF9Xk4gMVwhXCExX3t6X2kgPSAxfSBcZnJhY3tcdGF1XzEoeF9pLW1fMSleMn17Mn1cQmlnZykgClxlbmR7YWxpZ24qfSQkCgoKUGFyIHByb3ByacOpdMOpIGRlIGNvbmp1Z2Fpc29uIGRlIGxvaSwgb24gdm9pdCBhbG9ycyBxdWUgCiQkbV8xfFgsWixcdGF1LHAgXHNpbSBcbWF0aGNhbHtOfVxCaWdnKFxmcmFje1xmcmFje1xtdV8wfXtcc2lnbWFfMF4yfSArXHRhdV8xIFxzdW1fe2kgPSAxfV5OIDFcIVwhMV97el9pID0gMX14X2kgfXtcZnJhY3sxfXtcc2lnbWFfMF4yfStcdGF1XzEgXHN1bV97aSA9IDF9Xk4gMVwhXCExX3t6X2kgPSAxfSB9ICxcZnJhY3sxfXsgXGZyYWN7MX17XHNpZ21hXzBeMn0rXHRhdV8xIFxzdW1fe2kgPSAxfV5OIDFcIVwhMV97el9pID0gMX19IFxCaWdnKSAkJAoKRXQgZGUgbcOqbWUgcG91ciB0b3VzIGxlcyBhdXRyZXMgcGFyYW3DqHRyZXMuIENlbGEgaW5jaXRlIMOgIHJlcHJlbmRyZSB0b3V0IG5vdHJlIGNvZGUgcG91ciBjYWxjdWxlciBzZXVsZW1lbnQgY2VsYSDDoCBjaGFxdWUgw6l0YXBlLCBzYW5zIGRldm9pciByZWNhbGN1bGVyIGxhIHZyYWlzZW1ibGFuY2UuCgpWb2lsw6AgdW4gYm9uIGV4ZXJjaWNlIHBvdXIgdm90cmUgcHJvamV0ICEKCiMjIEV0IGF2ZWMgbGVzIHBhY2thZ2VzIAoKCgpgYGB7cn0Kc2V0d2QoIi9ob21lL2VzcGluYXNzZS9CdXJlYXUvVEFGL1NlYWZpbGUvQ291cnMvc3RhdCBiYXllc2llbm5lcy90cCBSIDIwMjAvb2xkX3NpdGUiKQpsaWJyYXJ5KHJqYWdzKQpOID0gMTAwMAp0cnVlID0gbGlzdChwID0gMC42LG0xID0gMCxtMiA9IDIsdGF1MSA9MSx0YXUyID0gNSkjI3BvdXIgbGEgc2ltdWxhdGlvbgp1ID0gcm1peCgxMDAwMCx0cnVlJHAsdHJ1ZSRtMSx0cnVlJG0yLHRydWUkdGF1MSx0cnVlJHRhdTIpCm0gPC0gamFncy5tb2RlbCgibW9kZWxfbWl4dHVyZS5idWciLCBkYXRhPWxpc3QoeCA9IHUsIE49TikpCgpmID0gY29kYS5zYW1wbGVzKG0sYygibSIsInRhdSIsInAiKSxuLml0ZXIgPSAxMDAwKQp0cmFqID0gZltbMV1dCnN1bW1hcnkoZikKcGFyKG1mcm93ID0gYyg1LDIpKQojIGZvciAoaSBpbiAxOjUpewojICAgcGxvdCgxOm5yb3codHJhaiksdHJhalssaV0sImwiLG1haW4gPSBjb2xuYW1lcyh0cmFqKVtpXSkKIyB9CnBsb3QoZikKYGBgCgoK